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Abstract 

We discuss a gauge-invariant prescription to take the mean-field approximation self- 
consistently in the PNJL model (Nambu-Jona-Lasinio model with the Polyakov 
loop) . We first address the problem of non- vanishing color density in normal quark 
matter, which is an artifact arising from gauge-fixed treatment of the Polyakov 
loop mean-fields. We then confirm that the gauge average incorporated in our pre- 
scription resolves this problem and ensures color neutrality. We point out that the 
proposed method has an advantage in computing the expectation value of any func- 
tion of the Polyakov loop matrix. We discuss the Casimir scaling as an immediate 
application of the method. 



Introduction The interplay between the QCD phase transitions of chiral 
restoration and color deconfinement at finite temperature and/or density has 
been attracting much interest recently. There are a lot of attempts to describe 
confinement-deconfinement physics by means of effective models in terms 
of the Polyakov loop [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]. 
One successful approach that can describe both the chiral and deconfinement 
transitions (crossovers) is the Polyakov-loop augmented Nambu-Jona-Lasinio 
(PNJL) model. This model accommodates self-consistent treatment for two 
approximate order parameters; the Polyakov loop L for deconfinement and 
the chiral condensate (ipip) for chiral restoration. Here the former works as an 
exact order parameter in the quenched limit (i.e. m q — > oo), while the latter 
is exact in the chiral limit (i.e. m q — > 0). Due to a particular form of coupling 
between L and (ipip) the PNJL model has a general tendency to make two 
crossovers in L and (ijjip) come close to each other [3]. Besides, it has turned 
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out that the bulk thermodynamics resulting from the model shows remarkable 
agreement with numerical data from the lattice QCD simulation [7]. 

It is known by now that pathological behavior arises from a simple mean- 
field ansatz for the Polyakov loop matrix [10,15,16,22]. In Ref. [23] one of the 
present authors found that a saddle-point approximation on the Polyakov loop 
matrix leads to unphysical non-zero color density even in the normal phase of 
quark matter (see also Refs. [24,25]). This is not a principle problem inherent 
to the PNJL model but rather a practical one associated with the mean-field 
approximation; we have to assume a certain gauge to make the color density 
definite, and at the same time, for the sake of the color density computation it 
is convenient to take a special gauge in which the Polyakov loop L is diagonal. 
These two gauge choices are, however, not necessarily compatible. In other 
words, the color chemical potential matrix and the Polyakov loop matrix are 
not commutable. 

The problem comes from the fact that we need to treat the Polyakov loop 
mean-field not as a traced quantity £ = -^-(trL) but as a matrix L((p 1 ,(p 2 ) = 
diag(e 1</31 , e 1</32 , e _1 ^ 1+ip2 - ) ) when we evaluate the color density. Then we face 
another undesirable situation. That is, in terms of ip± and y?2, it is hard to 
realize a difference between the Polyakov loop £ and the anti-Polyakov loop 
i = -^-(tiL^) which are both real numbers [11,26]. It is claimed in Ref. [18] 
that the fluctuation around the mean-field induces a difference between £ and 
£. We should note, however, that all these problems do not appear if we treat 
£ and £ as the relevant mean-fields, which works unless we consider color 
degrees of freedom such as color superconductivity [27]. In fact, in the color- 
superconducting phase, we cannot express the quasi-quark contribution to the 
thermodynamic potential solely in terms of £ and £, but it inevitably involves 
the matrix elements of L [10]. 

The present Letter aims to propose a resolution to circumvent these short- 
comings of the simple mean-field approximation. We would emphasize that 
our prescription not only improves the mean-field approximation but also en- 
compasses correct gauge dynamics from which the neutrality with respect to 
gauge charge is derived (i.e. the Gauss law). 



Model and mean-field approximation We here explain our model, the 
ingredients of which are the Polyakov-loop matrix model [5,9] and the NJL 
model [28] . The difference from the PNJL model lies in a mean- field evaluation 
for the Polyakov loop [1,3,6,29]. 

First, let us address the pure gluonic sector. We assume that the pure gluonic 
dynamics would be described by the nearest neighbor interaction of the traced 
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Polyakov loop as 



S e [L] = -N*e- a / T j:i(x)l*(x + h), (1) 

with / = ]^trL and I* = -^-trl/t. This form of the simplest matrix model [30] 
is to be postulated from the leading-order contribution in the strong coupling 
expansion, which specifies the T-dependent interaction strength with a being 
a model parameter [3] . 

The action (1) looks like a spin model. We then make use of the Weiss approx- 
imation with the neighboring spin sites treated as the mean-fields. Hence, the 
mean-field action is 

Srffafil = -N c J2[aRel(x) +i/3ImZ(x)], (2) 

X 

where a and (3 correspond to the Polyakov loop mean-fields. Finite (3 would be 
induced by C-odd terms at finite fi. We denote the Polyakov loop expectation 
values, hereafter, as £ = (l) m t and £ = (l*) m f. The expectation value (. . . ) mf 
refers to the average over the Polyakov loop matrix with the mean-field action. 

In the Weiss mean-field approximation the free energy is defined byPI 

= (S g [L\ - S m{ [L]) m{ -In J VLe- s ^ L \ (3) 

It is possible to find a closed analytic expression of / g (a, (3 = 0) at fi = 0, but 
we have to rely on numerical calculation to evaluate f g (a,f3) for fi ^ 0. 

We can fix the parameter a by requiring that the pure gluonic theory has 
a first-order phase transition of color deconfinement at T = 270 MeV when 
/i = 0. This condition results in 

a = 542 MeV. (4) 



Next, we shall consider how to add the contribution of dynamical quarks 
in the mean-field approximation. We simply add dynamical quarks using the 
quasi-quark approximation with the same Polyakov loop coupling as the PNJL 
model. In our notation ^rfi q (a"j, L) denotes the quark thermodynamic potential 
with the chiral condensates, a u , ad, and a s , giving the total mean- field free 
energy, 

f mf (<T i5 a, (3) = / g (a, (3) + (ftqfo; £)) mf - (5) 

1 Although the variational principle seems to break down due to the sign problem at 
fj, 0, the saddle-point of this mean-field free energy leads to a good approximation. 
See Ref. [11] for details. 
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Although f2 q is a complex function of L for /i 7^ 0, its expectation value as a 
function of a and (3 is real. This is because the imaginary part contributing 
to the thermodynamic potential is odd under C, i.e. L — > transformation. 
That is, 



1 r 

ft q (<7i,L)) = / dLe 



N c aRe I 



cos(iV c /3Im /)Re fi q - sin(iV c /3Im /)Im Q c 



(6) 



where z m { is the normalization given as z m f = J dL Q N ^ aRel cos(N c /3lm I), which 
is manifestly real. As for the Polyakov loop, we readily find 



— / dLe N < aRel 

z mf J 

— I dLe N < aRel 



cos(iV c /3Im/)Re/ - sin(iV c /3ImZ)ImZ 
cos(iV c /3Im l)Re I + sin(iV c /?Im Z)Im I 



(7) 



It is obvious from the above that i and £ are different by the presence of the 
imaginary (C-odd) part induced by (3 7^ at finite fi [11,31]. 



We now must specify the concrete form of Q q (<Ji,L). To this end, here, we 
shall augment the NJL model with the Polyakov loop coupling (i.e. the PNJL 
model). Then fl q (ai,L) take the following form; 

n q (a t , L) = g s {a 2 u + a 2 d + a 2 s ) + Ag D a u a d a s - 2N C £ J -^-^ £i (p) 



- 2T EE/ ^3 lndet ( 1 + LAe " (£ 



i X=±l ■ 



where the quasi-quark dispersion relations are Si{p) = yp 2 + M? with the 
constituent quark masses being M u = m u — 2gscr u — 2gr)0'd<Js, Md = ma — 
2g s a d - 2g D cr s a u , and M s = m s - 2g s a s - 2g D a u a d . We note that A = +1 
and —1 in the above are the quasi-quark and quasi- antiquark contributions, 
respectively. 

We take the same parameter set in the NJL part as in Ref. [28]; A = 631.4 MeV, 
m u = m d = 5.5 MeV, m s = 135.7 MeV, g s k 2 = 3.67, and g D A 5 = -9.29. Then 
&u = &d always holds due to isospin symmetry in the strong interaction. 

This model has one more parameter, that is the normalization of / g (a, 0). The 
right-hand side of Eq. (3) is proportional to the number of space points, N, 
and thus f g (a,{3) oc T • N/V which carries the mass dimension of the energy 
density. Here, N/V is a model parameter corresponding to b discussed in 
Ref. [21]. We can fix N/V by the condition that the chiral and deconfinement 
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crossovers take place near T = 200 MeV. In this way, we find 

N/V = 0.02A 3 . 



(10) 



Quark and color densities Once we determine the mean-fields {a, (3, <Ji] 
by solving the gap equations, df m f/da = df m f/dj3 = df m {/d<7i = 0, we can 
calculate various physical quantities. 



The quark number density, i.e. n q = —df m f/d/i, can be expressed as 



n 



1 ^ r d 3 p 



x tr 



J dL e ^( aRei +i/3Imi) 



1 



Lt e ( £ i(p)-^)/ T + 1 LefeW+^ + l 



Next we consider color densities. Since the phase A 4 of the Polyakov loop 
matrix, L = exp[L4 4 /T], could be regarded as the color chemical potential, 
the color density is then given by differentiating the integrand of / m f with 
respect to {—'\Af). This leads to 



_ J_V f ( dT p JVc(aRei+i/3Im0 

a ~ z^^J (2ttW 



2mf 

x tr 



1 



T — 



1 



Xte( e i(p)-M)/ T + 1 " a L e<f s *W+M T + 1 Ta 



'12) 



Here T a 's are the SU(iV c ) algebra in the fundamental representation. In de- 
riving this we have made use of the cyclicity in the trace. The group integra- 
tion in Eq. (12) is hard to perform in general. Usually we take the Polyakov 
gauge in which L is diagonal with two angle variables tpi and ip 2 , and then 
we can express dL as d(pidip 2 accompanied by the SU(iV c = 3) Haar mea- 
sure fj,(ipi,ip 2 ) = [sin(y?i - cp 2 ) + sin(2<y9i +<^ 2 ) + sin(y>i + 2(f 2 )} 2 / (3tt 2 ) . This 
procedure works straightforwardly for Eq. (11) but not for Eq. (12) because 
of the presence of T a . The color density is gauge dependent, however, so we 
should fix the gauge to define this quantity. If we take the Polyakov gauge as 
usual, then the color density for only the T 3 and T 8 components (belonging 
to the Cartan subalgebra of SU(3)) have non- vanishing integrands. We then 
can define the red, green, and blue quark densities as n r = \n q + ^n 3 + 



3 U 1 



~ l^ 3 + 273 n§ ' and Uh 



\n q 



"8, 



73 



n 8 . 



Approximation Here we introduce an approximation which greatly reduces 
the computational cost. That is, 



lndetfl + Le 



mf 



mf 



(13) 
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and the same for the antiquark part. With this approximation applied to the 
free energy expression, we can reduce the three-dimensional integral with re- 
spect to {</?i, <f2,p} to the one-dimensional p-integral with given £ and £ which 
result from the integral over {ipi, ip 2 } independently of p. We have numerically 
confirmed that this approximation works excellently well. 



Standard PNJL model treatment For comparison to the simple mean- 
field approximation used in literature [10] , we shall calculate the same physical 
quantities using the standard PNJL model, 

/mf(Ci, L) = V g i ue [l, I] 

/a d^p 
(2tx) 6 

2T E E | °^lndet(l+L A e-^( p )- A ^ T ), (14) 



i A=±l' 



where 

V gUe [l,l] = -6T|54e- a/T /r+ln[l-6/r-3(/0 2 + 4(/ 3 + P)]| (15) 

with a = 664 MeV, b = 0.026A 3 so that the transition (crossover) temperatures 
with and without dynamical quarks are 270 MeV and 200 MeV respectively, as 
explained previously. We will refer to this model as the standard PNJL model 
hereafter. 

The above expression for / m f(o"j,L) is given in terms of the gauge invariant 
mean-fields, £ and £, so that one can evaluate them without difficulty to find 
i > £ at non-zero fi. The serious problem arises when we are interested in 
quantities associated with color degrees of freedom. We cannot express the 
color density, n a , using £ and £ only, as seen from Eq. (12). That is also the 
case if the color-superconducting phase is considered in the PNJL model. 

Practically, in such a situation, one may well assume the mean-fields, ipi 
and </? 2 (or 3 and </>§), to characterize the Polyakov loop matrix as L = 
diag(e i¥,1 ,e i¥,2 ,e- i( ^ 1+ ^ 2) ) = exp[i(0 3 T 3 + 8 T 8 )/T]. This prescription is quite 
problematic, however, though adopted frequently. The traced Polyakov loops, 
£ and £, become complex with non-zero 3 and <j)$ in general. To avoid this 
artifact, one could assume 0g = 0, but such an assumption is not consistent 
with £ £ at finite /i. More seriously an unphysical color density is induced 
by the mean-fields, which is, of course, an artifact of this prescription. One 
may want to cancel the color density by introducing color chemical potentials 
[23,24], but as soon as one does so, another undesirable problem seems to come 
out immediately [32]. 
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Fig. 1. The chiral and deconfinement crossovers at /x = 0. The constituent quark 
masses are normalized by the vacuum values, M u q and M s q. The thin lines show 
the results from the standard PNJL model given in Eq. (14). 

We define the magnitude of the color density by n c = (ELi n a)^ 2 [23]. This 

quantity is invariant under gauge rotation. Thus, under the assumption 8 = 0, 

the color density magnitude is n c = -2=\n((f)s) — n(0)\, where we define 

v 3 



n[<P3 



d 3 p 



Re 



1 



1 



3 (ei(p)-M-i<fo)/T _|_ X e (£,(p)+M+i03)/T _j_ X 



(16) 



We remark that n c is no n- vanishing at finite 03. This is simply because the 
phase of the Polyakov loop matrix generally has the physical meaning of the 
color imaginary chemical potential, and so 03 ^ induces n c ^ 0. 



Order parameters and color densities We first show the order param- 
eters for chiral restoration and color deconfinement in Fig. 1 as a function of 
T at ix — 0. The thick curves represent the constituent quark masses and the 
Polyakov loop obtained from Eq. (5), while the thin curves are the results in 
the standard PNJL model with Eq. (14). We note that the thick and thin 
curves are very close, which means that our formulation would not spoil the 
nice feature established in the standard PNJL model. 

Now we shall move on to the finite density case. In Fig. 2 we display the 
physical quantities as a function of T at \i = 300 MeV. The left figure shows 
the constituent quark masses and the Polyakov loop in the same way as in 
Fig. 1. The thin lines are the results from the standard PNJL model again. 
The right figure shows the quark number density n q and the color densities 
{n r ,n g ,n b }, where n q = n r + n g + n& should be fulfilled. We see that the 
thick curves have significant difference from the thin curves resulting from the 
standard PNJL model. In the case of our prescription we have n r = n g = rib, 
meaning that n c = \n r + n g - 2n b |/ v / 3 = 0, while n c is non-vanishing in the 
standard PNJL model as depicted by the thin line with the label n c . It is 
also notable that, despite drastic difference in the color densities, the quark 
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Fig. 2. Left: Constituent quark masses, M u and M s , Polyakov loops, £ and £, as a 
function of T at \i = 300 MeV. The thin curves are the results from the standard 
PNJL model. Right: Color densities {n r ,n g ,rib} and the quark number density n q 
as a function of T. The net color density n c = \n r + n g — 2re&|/-\/3 is indicated by 
the thin line with the label n c in the case of the standard PNJL model. In our 
prescription n c is always zero. 

number density, n q , hardly changes; the thin line stays close to the thick line 
for n q . 



Casimir scaling at finite density As already noted, the mean-field ap- 
proximation discussed here enables us to compute not only the traced Polyakov 
loop in the fundamental representation but also the expectation value of arbi- 
trary functions of the Polyakov loop. We shall take a close look at the Polyakov 
loop in the higher representations as an immediate application. 

The Polyakov loop in the higher representations is of special interest with 
regard to the Casimir scaling hypothesis [33], which may provide a crucial key 
to understanding non-perturbative aspects of QCD such as confinement [30]. 
The Casimir scaling hypothesis claims that the color singlet potential between 
static color sources in the representation r is proportional to the Casimir 
invariant C^r). The statement is rather obvious in the perturbative regime, 
but it is quite non-trivial at large distances. From the theoretical perspective 
this hypothesis is verified up to two-loop order in the lattice perturbation 
theory both in pure gauge theory [34] and in QCD with massless dynamical 
quarks [35]. Beyond two-loop order the Casimir scaling can be violated, though 
the violation is tiny [34]. 

The scaling hypothesis is also tested numerically in the lattice simulation. 
In the SU(3) pure gauge theory at T = the hypothesis has been verified 
up to the string breaking distance [33]. The Casimir scaling hypothesis also 
brings a strong constraint on the Polyakov loop expectation value; the traced 
Polyakov loop in any representation r should satisfy the following scaling 
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Table 1 

Group theoretical factors in various representations; r is the representation labeled 
by its dimension D(r), (p,q) is the corresponding weight factor, t = p — q modulo 
3 is the triality, z = e l27T / 3 is an element of Z3, 62(7") is the quadratic Casimir 
invariant, and d r defines the ratio d r = C2(r)/C2(3). The dimension is given by 
D{r) = (p + l)(g + l)(p + q + 2)/2. In the triality zero representation (i.e. z = 1), 
the Polyakov loop is insensitive to center symmetry, and thus it does not serve as 
an order parameter of deconfinement. 

relation irrespective of the renormalization of the Polyakov loop [30]; 

ej dr « £ 3 . (17) 

Here £3 = I is the Polyakov loop in the fundamental representation, d r is the 
ratio of the quadratic Casimir invariant; d r = C2(r)/C2(3) = |C2(r). The 
relation (17) between the Polyakov loops in different representations actually 
provides a useful test for the hypothesis, and in Ref. [30] this test has been 
extensively performed with use of lattice QCD data both in pure gauge theory 
and in Nf = 2 QCD. It has been found that the scaling violation is visible 
only in the very vicinity of the first-order phase transition in pure gauge the- 
ory, while the deviation from the scaling is more evident in Nf = 2 QCD 
particularly below the crossover temperature. The scaling (17) is almost exact 
at high temperature where I is substantially large. 

It has been reported recently that the Casimir scaling of the Polyakov loop 
in the fundamental (3) and adjoint (8) representations is well realized in the 
PNJL model at \i — [36]. We here present the first systematic model study on 
the Casimir scaling at non-zero chemical potential. We compute the Polyakov 
loop in various representation from 3 to 27 as listed in Tab. 1. We can con- 
struct the Polyakov loop matrix in higher representations by the direct prod- 
ucts of L and V s in the fundamental representation using the Clebsch-Gordan 
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50 100 150 200 250 300 100 200 300 400 500 

T [MeV] |i [MeV] 

Fig. 3. Left: Scaled Polyakov loop l^ dr as a function of T for fx = and for 
H = 300 MeV in various representations with d r = C2(r)/C2(3) being the ratio 
of the quadratic Casimir invariant. Right: Scaled Polyakov loop as a function of /x 
for T = 100 MeV and for T = 150 MeV. It should be noted that the vertical axis is 
logarithmic. 

coefficients. In the left of Fig. 3 we show the various Polyakov loops as a func- 
tion of T at fi = 300 MeV as well as at \x = 0. We see that the scaling is good at 
high T in both cases. The scaling regime is reached faster in the \x = case; the 
violation of the Casimir scaling is small at the crossover around T = 200 MeV, 
while in the case of /j, = 300 MeV the deviations in l x J dr are significant at the 
crossover temperature T = 120 MeV. The presence of finite density tends to 
enhance the scaling violation. At the same time we should be careful about the 
interpretation; in Ref. [30] it has been shown that this kind of matrix-based 
model fails in reproducing the exact Casimir scaling £ x J dr = £3 observed on 
the lattice in the pure gauge theory. One comment which we should mention 
here is that the scaling violation at higher representations may be attributed 
to the fact that we limit ourselves to the simplest version of the matrix model 
in Eq. (1). It could be possible that the inclusion of the Polyakov loops in 
higher representations may diminish artificial violation of the Casimir scaling. 
For example we could consider, l 3 (x)l 6 (x + n), / 6 (x)/g(x + n), etc., in the model 
action, which are allowed by Z3 center symmetry in the pure gluonic sector. 
Moreover, as for the Polyakov loops in the triality zero representations such as 
8, 10, and 27, one may well add their arbitrary functions in the model action. 
It would be an interesting future problem to take account of the Polyakov loop 
in higher representations into the matrix model in such a way that the model 
preserves charge conjugation symmetry [31]. 

In the right of Fig. 3 we show the Polyakov loop l r as a function of [i for 
T = 100 MeV and for T = 150 MeV. At T = 100 MeV the first-order chiral 
transition occurs at fi — 315 MeV, while chiral restoration is smooth crossover 
with increasing /1 when T = 150 MeV. Again we notice the significant scaling 
violation, though the violation is exaggerated on the logarithmic plot. We 
observe in Fig. 3 that some £^ dr ^s cross each other as /1 increases. We see that 
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the change in 6 is milder than those in 3 and 15. This is reasonable, for the 
excitation with the triality z* (like 6) should be easier than that with z (like 
3 and 15) in a medium carrying the triality z at \i > 0. 



Conclusion We showed that the pathological problems associated with a 
simple mean-field approximation in the PNJL model can be resolved by the use 
of the Weiss mean-field approximation. We explicitly demonstrated that the 
color density is vanishing in the normal phase of quark matter. This vanishing 
color density is guaranteed by the integration with respect to the Polyakov 
loop, which can translate into the Gauss law resulting from the ^-integration. 
We also confirmed that £ > £ at finite \i is naturally realized. 

Our mean-field prescription allows us to compute the expectation value of any 
function of the Polyakov loop matrix, L and U , easily. As a demonstration, 
we computed the Polyakov loops in various representations from 3 to 27. We 
observed that the Casimir scaling is violated at finite /i more than at zero 
density, which is quite natural. More interestingly, the Polyakov loop with the 
triality z turned out to have decreasing behavior as a function of \i as long 
as is small. This means that the quark excitation bearing the same triality 
charge with the background medium is less favored. It would be an interesting 
future work to incorporate couplings between the Polyakov loops in higher 
representations into the matrix model and investigate the scaling violation in 
wider model space. 

An interesting extension of the present work would be the QCD phase struc- 
ture with the diquark condensation taken into account. It would be of par- 
ticular interest how the diquark condensates and the colored Polyakov loop 
matrix are entangled in a color-superconducting medium. The work along this 
line certainly deserves future investigations. 
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